Examples with lm() models.
But glm() (generalized linear models: Poisson, logistic, binomial, …) works similarly

summary table after fitting a model with lm()

1. Hypothesis testing
2. Regression
3. ANOVA
4. ANCOVA
5. Mixed models

1. Hypothesis testing
2. Regression
3. ANOVA
4. ANCOVA
5. Mixed models

• assumptions: random variable \(X\sim\text{normal}(\mu,\sigma)\)
• data \(x_1,\dots,x_n\)
• H0: Null hypothesis with random variable \(X\) - \(\mu=\text{E}(X)=0\)
• H1: Alternative hypothesis - e.g. \(\mu>0\)
• Test statistic \(T(X)\) transformation
  distribution of \(T\) under H0 is known - \(T(X)=\frac{\bar X-\mu}{\sigma}=\frac{1}{\sigma}\frac{1}{n}\sum_{i=1}^n{x_i}\sim\text{normal}(0,1)\)

• define type I error \(\alpha\) (e.g. \(\alpha=0.05\)):
  chance of rejecting H0 (\(\mu=0\)) although it is true
• check if \(T(\text{data})\) is outside \((1-\alpha)\) quantile
  ➔ reject H0!
• P-value: \(P(T>T(\text{data}))\)

\(X_i\) have normal distribution

\(Z=\sum_i X_i\) sum has normal distribution

e.g. Z-test for mean (standard deviation \(\sigma\) is known)

\(X_i\) have standard normal distribution

\(Z=\sum_iX_i^2\) sum of squares has \(\chi^2_n\) distribution

e.g. chisq-test for independence

\(X\sim\chi^2_n\), \(Y\sim\chi^2_m\)

\(Z=\displaystyle\frac{X/n}{Y/m}\) ratio has \(F_{n,m}\) distribution

e.g. F-test for regression, ANOVA, for comparing variances

\(X\) normal, \(Y\sim\chi_n^2\)

\(Z=\displaystyle\frac{X}{\sqrt{Y/n}}\) has \(t_n\) distribution

e.g. t-test for regression coefficients, 2-sample t-test for means

library(palmerpenguins)
data=as.data.frame(penguins)
str(data)

## 'data.frame':    344 obs. of  8 variables:
##  $ species          : Factor w/ 3 levels "Adelie","Chinstrap",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ island           : Factor w/ 3 levels "Biscoe","Dream",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ bill_length_mm   : num  39.1 39.5 40.3 NA 36.7 39.3 38.9 39.2 34.1 42 ...
##  $ bill_depth_mm    : num  18.7 17.4 18 NA 19.3 20.6 17.8 19.6 18.1 20.2 ...
##  $ flipper_length_mm: int  181 186 195 NA 193 190 181 195 193 190 ...
##  $ body_mass_g      : int  3750 3800 3250 NA 3450 3650 3625 4675 3475 4250 ...
##  $ sex              : Factor w/ 2 levels "female","male": 2 1 1 NA 1 2 1 2 NA NA ...
##  $ year             : int  2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 ...

2 sample t-test

Data in 2 groups: \(x_1\dots x_n\)
                               \(y_1\dots y_m\)

Assumption: \(X\sim\text{normal}(\mu_X,\sigma_X)\)
                        \(Y\sim\text{normal}(\mu_Y,\sigma_Y)\)

H0: \(\mu_X=\mu_Y\) (means are identical)

H1: \(\mu_X\neq\mu_Y\)

under H0: \(T=( \bar X - \bar Y ) / \sqrt{\frac{s^2_X}{n}+\frac{s^2_Y}{m}}\) has t-distribution

x = subset(data, species=="Chinstrap")$flipper_length_mm
y = subset(data, species=="Adelie")$flipper_length_mm
boxplot(x,y)

t.test(x,y)

## 
##  Welch Two Sample t-test
## 
## data:  x and y
## t = 5.7804, df = 119.68, p-value = 6.049e-08
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  3.859244 7.880530
## sample estimates:
## mean of x mean of y 
##  195.8235  189.9536

➔ reject H0 (\(P<0.05\)) that the means are identical

## t = 5.7804, df = 119.68, p-value = 6.049e-08

curve(dt(x, df=119.68), from=0, to=6, lwd=3, col=4)
abline(v=qt(p=0.95, df=119.68), lty=2  )
points(5.7804, 0, col=2, pch=16, cex=1.5)

1. Hypothesis testing
2. Regression
3. ANOVA
4. ANCOVA
5. Mixed models

summary table after fitting a model with lm()

\(SST=\sum_i({y}_i-\bar{y})^2\)     total
\(SSE=\sum_i(\widehat{y}_i-\bar{y})^2\)    explained
\(SSR=\sum_i({y}_i-\widehat{y}_i)^2\)   residual

\(SST=SSE+SSR\)

variation explained by the model

\(R^2=\frac{SSE}{SST}=1-\frac{SSR}{SST}\)

• two nested models, e.g.    m1: y~x1*x2
                                               m2: y~x1+x2
• m1 better than m2?
  residuals: \(SSR_1<SSR_2\)?
• H0: \(SSR_1=SSR_2\)
• test statistic \(\frac{(SSR_2-SSR_1)/(k_1-k_2)}{SSR_1/(n-k_2-1)}\) has \(F_{k_1-k_2,n-k_2-1}\) distr.
  (ratio of two \(\chi^2\) distributions)
• If \(P<0.05\) ➔ reject H0
• H1: m1 has better fit

• test your regression model     y~1+x1+...+xk
  against intercept only model  y~1
• H0: regression model is only as good as intercept model
         same as \(SSR=SST\)
         same as \(R^2=0\)
         same as \(\beta_1=\ldots=\beta_k=0\)
• test statistic \(\frac{SSE/k}{SSR/(n-k-1)}\) has \(F_{k,n-k-1}\) distribution
• If \(P<0.05\) ➔ reject H0
• H1: at least one effect \(\beta_j\neq0\)
  regression model better than intercept only

## Multiple R-squared: 0.4421, Adjusted R-squared: 0.416
## F-statistic: 16.91 on 3 and 64 DF, p-value: 3.389e-08

curve(df(x,3,64), from=0, to=17, col=4, lwd=2)
abline(v=qf(p=0.95, df1=3, df2=64), lty=2)
points(16.91, 0, col=2, pch=16, cex=1.5)

• test full regression model   y~x1+x2+x3
  against reduced model       y~ x2+x3
• H0: full model is only as good as reduced model
         same as \(SSR_f=SSR_r\)
         same as \(\beta_1=0\)
  
• F test for \(SSR_f=SSR_r\)    ➔ \(F_{1,n-k-1}\)
  t test for \(\beta_1 / se(\beta_1)=0\)         ➔ \(t_{n-k-1}\)       (identical)
• If \(P<0.05\) ➔ reject H0
• H1: \(\beta_1\neq0\)
• ATTN: if x1 involved in interaction x1:x2,
              test for main effect \(\beta_1\) is meaningless

## Estimate Std. Error t value Pr(>|t|)
## body_mass_g:bill_length_mm 3.371e-04 5.647e-04 0.597 0.5527

curve(dt(x,64), from=-4, to=4, col=4, lwd=2)
abline(v=qt(p=0.975, df=64), lty=2) # two-sided: alpha/2 = 2.5%
points(0.597, 0, col=2, pch=16, cex=1.5)

multiple regression y~1+x1+...+xk

• H0: model 1 and model 2 perform equally (\(SSR_1=SSR_2\))
• If \(P<0.05\) reject H0
• H1: model 2 performs better than model 1 (\(SSR_1>SSR_2\))

multiple regression y~1+x1+...+xk

• stepwise inclusion of predictors
  
• interactions x1:x2 treated as individual predictor
  
• order matters for table: y~x1+x2 or y~x2+x1

anova(mod1)

## Analysis of Variance Table
## 
## Response: flipper_length_mm
##                            Df  Sum Sq Mean Sq F value    Pr(>F)    
## body_mass_g                 1 1402.68 1402.68 47.2190 3.124e-09 ***
## bill_length_mm              1   93.44   93.44  3.1457   0.08089 .  
## body_mass_g:bill_length_mm  1   10.58   10.58  0.3563   0.55267    
## Residuals                  64 1901.17   29.71                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• 1 + body_mass_g better than 1
• 1 + body_mass_g + bill_length_mm
  not better than 1 + body_mass_g

summary.aov(mod1)

##                            Df Sum Sq Mean Sq F value   Pr(>F)    
## body_mass_g                 1 1402.7  1402.7  47.219 3.12e-09 ***
## bill_length_mm              1   93.4    93.4   3.146   0.0809 .  
## body_mass_g:bill_length_mm  1   10.6    10.6   0.356   0.5527    
## Residuals                  64 1901.2    29.7                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• same output

• Careful with anova(mod1)
• Type I anova() table usually not recommended in regression with multiple predictors
• rather use summary()
• anova(mod1, mod2) F-tests for competing models
• for ANOVA with categorical predictors
  Type II ANOVA often makes more sense

\(AIC\)   \(=\)   \(-2\log L\)   \(+\)   \(2k\)           lower=better
                  goodness-of-fit          penalty

• compare models, lower AIC = better
• tradeoff between goodness-of-fit and model complexity \(k\)
• models don’t need to be nested
• pure value doesn’t say anything (as opposed to \(R^2\)),
  only model difference \(AIC_1-AIC_2\)
• only compare models on exactly same data,
  watch out for NA in data

\(BIC\)   \(=\)   \(-2\log L\)   \(+\)   \(k\log n\)           lower=better
                  goodness-of-fit              penalty

• \(\log L=\sum_{i=1}^n \log lik(y_i)\) grows with data size \(n\)
• choose a penalty term that also grows with \(n\)
• penalizes model complexity \(k\) more strongly than \(AIC\)
• tends to select for more parsimonious model than \(AIC\)
• note: BIC does not require Bayesian statistics

# mod1 = lm(flipper_length_mm ~ body_mass_g * bill_length_mm, data)
# mod2 = lm(flipper_length_mm ~ body_mass_g , data)
AIC(mod1, mod2)

##      df      AIC
## mod1  5 429.4645
## mod2  4 427.8421

very small difference, still select simpler model

BIC(mod1, mod2)

##      df      BIC
## mod1  5 440.5621
## mod2  4 436.7201

BIC difference larger than AIC difference

• anova(mod1, mod2): F-test for comparing nested models
• summary.lm(mod1): F-test for overall model fit and \(R^2\)
• summary.lm(mod1): t-tests for single parameters
• summary.aov(mod1): table for stepwise F-tests
  ➔ not recommended for regression
• AIC(mod1, mod2) and BIC(mod1, mod2)

1. Hypothesis testing
2. Regression
3. ANOVA
4. ANCOVA
5. Mixed models

library(palmerpenguins)
data = as.data.frame(penguins)
boxplot(flipper_length_mm ~ species, data=data)

library(ggplot2)
ggplot(data=data, aes(y=flipper_length_mm, x=species)) + 
  geom_boxplot()

response y

one grouping factor A with several levels A1, A2, A3, …

test for difference between group means

H0: all group means identical, \(\bar{y}_1=\bar{y}_2=\bar{y}_3=\ldots\)
➔ identical to overall mean \(\bar{y}\)

F-test / ANOVA for y~1 vs y~A

\(P<0.05\) reject H0

➔ H1: not all group means identical

mod1 = lm(flipper_length_mm ~ species, data=data)
summary.aov(mod1)

##              Df Sum Sq Mean Sq F value Pr(>F)    
## species       2  52473   26237   594.8 <2e-16 ***
## Residuals   339  14953      44                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 2 observations deleted due to missingness

"y~species" is better than "y~1"

But how to get group-level means (predictions)?

Means using model.tables(). Requires aov() object.

aov1 = aov(flipper_length_mm ~ species, data=data)
model.tables(aov1, "means")

## Tables of means
## Grand mean
##          
## 200.9152 
## 
##  species 
##     Adelie Chinstrap Gentoo
##        190     195.8  217.2
## rep    151      68.0  123.0

Contrasts using TukeyHSD(). Requires aov() object.

# aov1 = aov(flipper_length_mm ~ species, data=data)
TukeyHSD(aov1)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = flipper_length_mm ~ species, data = data)
## 
## $species
##                       diff       lwr       upr p adj
## Chinstrap-Adelie  5.869887  3.586583  8.153191     0
## Gentoo-Adelie    27.233349 25.334376 29.132323     0
## Gentoo-Chinstrap 21.363462 19.000841 23.726084     0

TukeyHSD corrects P-values for multiple comparisons.

Plot means and data using sjPlot.

library("sjPlot")
plot_model(mod1, type="pred", show.data = TRUE, jitter=0.5)

emmeans() can predict means and also contrasts,
works with lm() object.

library(emmeans)
emmeans(mod1, ~species)

##  species   emmean    SE  df lower.CL upper.CL
##  Adelie       190 0.540 339      189      191
##  Chinstrap    196 0.805 339      194      197
##  Gentoo       217 0.599 339      216      218
## 
## Confidence level used: 0.95

Contrasts differences using emmeans().

pairs(emmeans(mod1, ~species))

##  contrast           estimate    SE  df t.ratio p.value
##  Adelie - Chinstrap    -5.87 0.970 339  -6.052  <.0001
##  Adelie - Gentoo      -27.23 0.807 339 -33.760  <.0001
##  Chinstrap - Gentoo   -21.36 1.004 339 -21.286  <.0001
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

Plot means using emmeans().

plot(emmeans(mod1, ~species), horizontal=FALSE)

• ANOVA is just a linear model
• F-test for difference in group-means
• y~A vs y~1
• lm() and aov() use dummy-coding.
  estimated parameters are not group-means
• predicted group-means (emmeans) are model-dependent!
  important in two-way ANOVA and GLM
• post-hoc tests / contrast:
  test for differences in certain groups (t-tests)
  P-values are corrected for multiple comparisons (e.g. Tukey)
• with a GLM, use anova( , test="Chi") for a Chi-squared likelihood-ratio-test (LRT) instead

1. Hypothesis testing
2. Regression
3. ANOVA
4. ANCOVA
5. Mixed models

data = subset(data, !is.na(sex))
ggplot(data = data, aes(y=flipper_length_mm, x=species, colour=sex)) + 
  geom_boxplot()

mod1 = lm(flipper_length_mm ~ species + sex, data=data)

estimated group-means:

• effect of sex independent of species
• effect of species independent of sex
• Df < level combinations

mod2 = lm(flipper_length_mm ~ species * sex, data=data)

estimated group-means:

• effect of sex changes with species
• effect of species changes with sex
• Df = level combinations
• fits independent means for all level combinations

two-way ANOVA w. interaction y~A*B, or y~1+A+B+A:B

• includes predictors stepwise
• tests against previous model
• order matters for the table

# mod2 = lm(flipper_length_mm ~ species * sex, data=data)
summary.aov(mod2) # same as anova(mod2)

##              Df Sum Sq Mean Sq F value  Pr(>F)    
## species       2  50526   25263 789.912 < 2e-16 ***
## sex           1   3906    3906 122.119 < 2e-16 ***
## species:sex   2    329     165   5.144 0.00631 ** 
## Residuals   327  10458      32                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

two-way ANOVA w. interaction y~A*B, or y~1+A+B+A:B

• independently removes predictors completely
  (with its higher-order interactions)
• tests against full model
• order does not matter for the table

# mod2 = lm(flipper_length_mm ~ species * sex, data=data)
library(car)
Anova(mod2) # same as default: Anova(mod2, type=2)

## Anova Table (Type II tests)
## 
## Response: flipper_length_mm
##             Sum Sq  Df  F value    Pr(>F)    
## species      50185   2 784.5829 < 2.2e-16 ***
## sex           3906   1 122.1189 < 2.2e-16 ***
## species:sex    329   2   5.1442  0.006314 ** 
## Residuals    10458 327                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

predicted means

# mod2 = lm(flipper_length_mm ~ species * sex, data=data)
library(emmeans)
emmeans(mod2, ~species*sex)

##  species   sex    emmean    SE  df lower.CL upper.CL
##  Adelie    female    188 0.662 327      186      189
##  Chinstrap female    192 0.970 327      190      194
##  Gentoo    female    213 0.743 327      211      214
##  Adelie    male      192 0.662 327      191      194
##  Chinstrap male      200 0.970 327      198      202
##  Gentoo    male      222 0.724 327      220      223
## 
## Confidence level used: 0.95

# emmeans(mod2, ~species|sex) # same, just ordered
# emmeans(mod2, ~sex|species) # same, just ordered the other way

in-group contrasts only

pairs( emmeans(mod2, ~species|sex) )

## sex = female:
##  contrast           estimate    SE  df t.ratio p.value
##  Adelie - Chinstrap    -3.94 1.174 327  -3.356  0.0025
##  Adelie - Gentoo      -24.91 0.995 327 -25.044  <.0001
##  Chinstrap - Gentoo   -20.97 1.221 327 -17.169  <.0001
## 
## sex = male:
##  contrast           estimate    SE  df t.ratio p.value
##  Adelie - Chinstrap    -7.50 1.174 327  -6.388  <.0001
##  Adelie - Gentoo      -29.13 0.981 327 -29.694  <.0001
##  Chinstrap - Gentoo   -21.63 1.210 327 -17.870  <.0001
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

in-group contrasts only

pairs( emmeans(mod2, ~sex|species) )

## species = Adelie:
##  contrast      estimate    SE  df t.ratio p.value
##  female - male    -4.62 0.936 327  -4.932  <.0001
## 
## species = Chinstrap:
##  contrast      estimate    SE  df t.ratio p.value
##  female - male    -8.18 1.372 327  -5.961  <.0001
## 
## species = Gentoo:
##  contrast      estimate    SE  df t.ratio p.value
##  female - male    -8.83 1.037 327  -8.518  <.0001

marginal contrasts

pairs( emmeans(mod2, ~sex) )

##  contrast      estimate    SE  df t.ratio p.value
##  female - male    -7.21 0.653 327 -11.046  <.0001
## 
## Results are averaged over the levels of: species

pairs( emmeans(mod2, ~species) )

##  contrast           estimate    SE  df t.ratio p.value
##  Adelie - Chinstrap    -5.72 0.830 327  -6.890  <.0001
##  Adelie - Gentoo      -27.02 0.699 327 -38.681  <.0001
##  Chinstrap - Gentoo   -21.30 0.860 327 -24.774  <.0001
## 
## Results are averaged over the levels of: sex 
## P value adjustment: tukey method for comparing a family of 3 estimates

marginal contrasts

mod2 |> emmeans(~sex) |> pairs()

##  contrast      estimate    SE  df t.ratio p.value
##  female - male    -7.21 0.653 327 -11.046  <.0001
## 
## Results are averaged over the levels of: species

mod2 |> emmeans(~species) |> pairs()

##  contrast           estimate    SE  df t.ratio p.value
##  Adelie - Chinstrap    -5.72 0.830 327  -6.890  <.0001
##  Adelie - Gentoo      -27.02 0.699 327 -38.681  <.0001
##  Chinstrap - Gentoo   -21.30 0.860 327 -24.774  <.0001
## 
## Results are averaged over the levels of: sex 
## P value adjustment: tukey method for comparing a family of 3 estimates

emmip(mod2, ~species|sex, CIs=TRUE) +
  geom_jitter( aes(x=species, y=flipper_length_mm), 
             data=data, width=0.1, alpha=0.2, color="red")

emmip(mod2, ~sex|species, CIs=TRUE) +
  geom_jitter( aes(x=sex, y=flipper_length_mm), 
             data=data, width=0.1, alpha=0.2, color="red")

• test if means change with factor levels
• y~A+B vs y~A*B
• Type 2 Anova() table for multiple predictors or interactions
• model-based group means with emmeans
  e.g. y~A+B and y~A*B have different means
• post-hoc tests / contrast (e.g. Tukey)
• with a GLM, Anova() automatically chooses a
  Chi-squared likelihood-ratio-test (LRT) instead

1. Hypothesis testing
2. Regression
3. ANOVA
4. ANCOVA
5. Mixed models

ggplot( aes(x=body_mass_g, y=flipper_length_mm, color=species), data=data) + 
  geom_point()

• individual intercepts per group (dummy-coded)
• slope identical across groups
• test for difference in group-means while correcting for covariate effect: y~A+x vs y~x

• individual intercepts per group (dummy-coded)
• individual slopes per group (dummy-coded)
• test for difference in slopes: y~A*x vs y~A+x

# mod2 = lm(flipper_length_mm ~ mass_z * species, data=data)
Anova(mod2) 

## Anova Table (Type II tests)
## 
## Response: flipper_length_mm
##                Sum Sq  Df F value  Pr(>F)    
## mass_z         5080.0   1 177.319 < 2e-16 ***
## species        5903.2   2 103.028 < 2e-16 ***
## mass_z:species  244.6   2   4.269 0.01478 *  
## Residuals      9368.2 327                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

estimated group-level slopes

# mod2 = lm(flipper_length_mm ~ mass_z * species, data=data)
library("emmeans")
emtrends(mod2, ~species, var="mass_z")

##  species   mass_z.trend    SE  df lower.CL upper.CL
##  Adelie            5.32 0.780 327     3.79     6.86
##  Chinstrap         9.59 1.370 327     6.89    12.28
##  Gentoo            7.52 0.791 327     5.97     9.08
## 
## Confidence level used: 0.95

estimated contrasts for group-level slopes

# mod2 = lm(flipper_length_mm ~ mass_z * species, data=data)
pairs( emtrends(mod2, ~species, var="mass_z") )

##  contrast           estimate   SE  df t.ratio p.value
##  Adelie - Chinstrap    -4.26 1.58 327  -2.704  0.0197
##  Adelie - Gentoo       -2.20 1.11 327  -1.978  0.1192
##  Chinstrap - Gentoo     2.06 1.58 327   1.305  0.3933
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

1. Hypothesis testing
2. Regression
3. ANOVA
4. ANCOVA
5. Mixed models

again, continuous predictor(s) x, categorical predictor A

• fixed effects only model:

     y~x+A separate intercepts model (ANCOVA)
     y~x*A separate intercepts & slopes model (ANCOVA)

• mixed effects model:

     y~x+(1|A) random intercepts model (LMM)
     y~x+(1+x|A) random intercepts & slopes model (LMM)

➔ How to find best model (fixed and random effects structure) with AIC / BIC / F-tests etc?

library(lme4)
lmer( y ~ x + (1+x|A), data=data, REML=TRUE )

Zuur et al.: Mixed effects models (2009)

1. Start with maximal fixed effects structure
   
2. Find optimal random effects structure using REML=TRUE
   
3. Find optimal fixed effects structure using REML=FALSE
   
4. Fit best model again using REML=TRUE

ANCOVA model from last section, add sex as random variable

lmmax1 = lmer( flipper_length_mm ~ mass_z * species + 
                 (1+mass_z|sex), data=data )

Fixed effects have maximal structure.

• Theory in Biodiversity Science group (TiBS) provides consultation in theory, maths and statistics.
• each Wednesday after the seminar
• write me an email with a brief description of your problem
  benjamin.rosenbaum@idiv.de